Summary work

Setting up the environment

# Packages
library(ggplot2)
library(dplyr)
library(plotly)
library(tidyr)
library(RColorBrewer)
library(limma)
library(gplots)
library(edgeR)
library(devtools)
library(ggbiplot)
library(DESeq)
library(geneplotter)
library(Rtsne)
set.seed(1234) # makes everything reproductible 
meta_data <- read.table("../Data/raw_data/meta_data.txt", header = TRUE)
dataInitial <- read.table("../Data/raw_data/expn_matrix.txt", header = TRUE, row.names = 1)
spike_data <- read.table("../Data/raw_data/expn_matrix_spike.txt", header = TRUE, row.names = 1)

Normalization

Kevin’s normalization

Here the data is normalized based on the spike_data we have. Although this seems to be correct, there could be some issues: See post.

spike_data <- rbind(spike_data, "Totals" = colSums(spike_data)) #finding sum of column and adding it as new row
scaled_totals <- spike_data["Totals", ]
average_reads <- rowMeans(scaled_totals)
scaled_totals <- rbind(scaled_totals, "Scale" = scaled_totals[1, ] / average_reads) #dividing total spike values by average spike values
dim_names <- list(row.names(dataInitial))       #setting dim_names for upcoming matrix
dim_names[2] <- list(colnames((dataInitial)))
normalized_data <- matrix(nrow = nrow(dataInitial), ncol = ncol(dataInitial), dimnames = dim_names) #pre-making a matrix
for(i in 1: ncol(dataInitial)){   #for loop which scales data and saves it as normalized_data
  normalized_data[ , i] <- dataInitial[ , i] / scaled_totals[2, i]
}
normalized_data <- as.data.frame(normalized_data)
#normalizing for cell size (cellular miRNA counts)
total_data <- rbind(spike_data, "Total Spike" = colSums(spike_data), "Total Reads" = colSums(dataInitial))
total_cellular <- total_data["Total Reads", ]
average_reads <- rowMeans(total_cellular)
cell_scales <- rbind(total_cellular, "Scale" = total_cellular[1, ] / average_reads) #dividing total spike values by average spike values
dim_names <- list(row.names(dataInitial))       #setting dim_names for upcoming matrix
dim_names[2] <- list(colnames((dataInitial)))
normalized_data2 <- matrix(nrow = nrow(dataInitial), ncol = ncol(dataInitial), dimnames = dim_names) #pre-making a matrix
for(i in 1: ncol(dataInitial)){   #for loop which scales data and saves it as normalized_data
  normalized_data2[ , i] <- normalized_data[ , i] / scaled_totals[2, i]
}
normalized_data2 <- as.data.frame(normalized_data2)
# deleting dead cells
alive_indexes <- list() #initializing lists
k <- 1  #initializing counter
for(i in 1:ncol(total_data)){
  if(total_data["Total Reads", i] > 20000){
    alive_indexes[k] <- i
    k <- k + 1
  }
}
alive_data <- normalized_data2[ , as.numeric(alive_indexes)]
alive_meta_data <- meta_data[as.numeric(alive_indexes), ]
#Saving completely normalized data as n_data and n_meta_data
n_data <- alive_data
n_meta_data <- alive_meta_data
#makes group and design matrix
data_day <- as.character(n_meta_data$Population)
data_day <- recode(data_day, "HL60D0" = "Day 0", "HL60D1" = "Day 1", "HL60D3" = "Day 3", "HL60D7" = "Day 7")
data_day <-factor(data_day)
design <- model.matrix(~0 + n_meta_data$Population)
#deletes "n_meta_data$Population" from name of columns
colnames(design) <- gsub("n_meta_data\\$Population", "", colnames(design))
#lets call intercept HL60D0
#colnames(design)[1] <- "HL60D0"
# Keeps genes without all zeros miRNA
dataNZ <- n_data[which(rowSums(n_data) > 0),] 
#making the contrast matrix
# let's suppose for now that we want to compare each group with each other ?????????????????????????
contrastMatrix <- makeContrasts(HL60D1-HL60D0,
                                HL60D3-HL60D0,
                                HL60D7-HL60D0,
                                HL60D1-HL60D3,
                                HL60D1-HL60D7,
                                HL60D3-HL60D7,
                                levels = c("HL60D0","HL60D1","HL60D3","HL60D7"))
#will not be needed: don't keep big data
remove(normalized_data)
remove(normalized_data2)
remove(alive_data)

Please note that with this normalization, we are only using cells ????????????????(alive cells?)

General normalization

To be sure that the normalization is good I will continue the analysis with a more general normalization. The following is based on this page. By general I mean not specific to pikes.

data_dayDESeq <- as.character(meta_data$Population)
data_dayDESeq <- recode(data_dayDESeq, "HL60D0" = "Day 0", "HL60D1" = "Day 1", "HL60D3" = "Day 3", "HL60D7" = "Day 7")
data_dayDESeq <-factor(data_dayDESeq)
deSeqDat <- newCountDataSet(dataInitial, data_dayDESeq)
# Note: actually it's not a real normalisation. Rather computing size factors depending on ratio of medians 
# if all size factors are roughly equal to one, the libraries have been sequenced equally deeply.
deSeqDat <- estimateSizeFactors(deSeqDat)
head(sizeFactors(deSeqDat))

HL60D0.C01_S1 HL60D0.C02_S2 HL60D0.C03_S3 HL60D0.C04_S4 HL60D0.C05_S5 HL60D0.C06_S6 
    0.5254052     1.1318130     0.9646220     0.3233751     1.3896389     0.9596427 

idx.nz <- apply(counts(deSeqDat), 1, function(x) { all(x > 0)})
nNZsamples <- sum(idx.nz)
#will not be needed: don't keep big data
remove(dataInitial)

We see that the number of non zero samples in all genes is very low compared to traditional RNA-seq: 5 nomrally in thousands. This may be a good reason to stick with the normalization with skie RNA.

#plotting the estimated dispersions against the mean normalized counts
deSeqDat <- estimateDispersions(deSeqDat)
plotDispEsts(deSeqDat)

multidensity( counts(deSeqDat, normalized = T),xlab="mean counts", xlim=c(0, 1000))

 multiecdf( counts(deSeqDat, normalized = T),xlab="mean counts", xlim=c(0, 1000))

The two charts above clearly shows that the second normalisation isn’t convincing. Indeed, the second chart assesses whether the normalization has worked, and the densities should overlapp since most of the genes are heavily affected by the experimental conditions. Note: The strange density chart could also be due to the fact that miRNA are very rarely expressed in every sample (as we have seen before).

Plots

PCA

PCA no standardization

First we use a log transformation to make the data approximatively follow the homoscedasticity assumption.

#Log transform
log_normalized_data <- log(n_data)
log_normalized_data[log_normalized_data == "-Inf"] <- 0
data_pca <- prcomp(t(log_normalized_data))
g <- ggbiplot(data_pca, 
              scale = 1, 
              obs.scale = 1, 
              varname.abbrev = FALSE,
              var.axes = FALSE,
              pc.biplot =TRUE,
              circle = TRUE, 
              groups = data_day, 
              ellipse= TRUE) +
  ggtitle("PCA without normilizing") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
g 

#will not be needed: don't keep big data
#remove(g)
remove(data_pca)
remove(log_normalized_data)

Note: The PCA graph is interesting as it shows that after each day the clusters seem to go down in PC2 and extend their variance in PC1.

PCA standardized

Knowing the importance of having standardized data for PCA, we could think of standarding ours. This in general is a good idea, as veraibles are often not on the same scale and thus cannot be compared directly. Here, on revanche we have dimensions that are comparable as they each represent the number of miRNA in each cell. It is therefore a good idea not to standrdize, as it will unnecessarily decrease the variance explained. More information on this can be found on the following forums: (here)[http://stats.stackexchange.com/questions/105592/not-normalizing-data-before-pca-gives-better-explained-variance-ratio] or (here)[http://stats.stackexchange.com/questions/69157/why-do-we-need-to-normalize-data-before-analysis]. Out of curiosity we will still plot the standardized data:

#Log transform
log_normalized_dataNZ <- log(dataNZ)
log_normalized_dataNZ[log_normalized_dataNZ == "-Inf"] <- 0
data_pcaNZ <- prcomp(t(log_normalized_dataNZ), center = TRUE, scale. = TRUE)
g <- ggbiplot(data_pcaNZ, 
              scale = 1, 
              obs.scale = 1, 
              varname.abbrev = FALSE,
              var.axes = FALSE,
              pc.biplot =TRUE,
              circle = TRUE, 
              groups = data_day, 
              ellipse= TRUE) +
  ggtitle("PCA non zero genes") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
g 

#will not be needed: don't keep big data
#remove(g)
remove(data_pcaNZ)
remove(log_normalized_dataNZ)

From this plot we see that as we’ve supposed, the variance explained is lower than in the first one. We thus conclude that we should stick with the first plot as we can (because on the same scale).

PCA top 50

Let’s now plot the top 50 genes

#arranging data based on total expression
indexTopFifty <- sort(rowSums(n_data), index=T, decreasing=TRUE)$ix[1:50]
topFifty <- n_data[indexTopFifty,]
# note: small difference with kenny probably come from the fact that he did log before and then took top fifty 
#Log transform
log_topFifty <- log(topFifty)
log_topFifty[log_topFifty == "-Inf"] <- 0
#PCA
data_topFifty <- prcomp(t(log_topFifty))
g <- ggbiplot(data_topFifty, 
              scale = 1, 
              obs.scale = 1, 
              varname.abbrev = FALSE,
              var.axes = FALSE,
              pc.biplot =TRUE,
              circle = TRUE, 
              groups = data_day, 
              ellipse= TRUE) +
  ggtitle("PCA top 50 genes") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
g

#will not be needed: don't keep big data
#remove(g)
remove(data_topFifty)
remove(log_topFifty)
remove(topFifty)

Correlation heatmap

# Note: both plots of dataNZ and full data gives same correlation heatmap
normalized_data_cor <- cor(dataNZ)  #calculating sample-sample correlations
scaled_normalized_data_cor <- scale(normalized_data_cor)  #scaling for better contrast
limit <- 2
scaled_normalized_data_cor[which(scaled_normalized_data_cor > limit)] <- limit    #trimming for better contrast
scaled_normalized_data_cor[which(scaled_normalized_data_cor < -limit)] <- -limit
mypal <- rev(colorRampPalette(brewer.pal(11,"RdBu"))(100)) #setting colour palette
alphabet <- c("#b2df8a","#33a02c","#ff7f00","#fdbf6f","#b15928","#1f78b4","#6a3d9a","#e31a1c","#a6cee3","#cab2d6","#ffff99","#fb9a99")
pops <- alphabet[as.factor(data_day)] #setting population colour palette
heatmapPlot <- heatmap.2(scaled_normalized_data_cor, Rowv = NA, Colv = NA, symm = T, #making a heatmap
          trace = "none", dendrogram = "none", col = mypal,
          labCol = NA, labRow = NA,
          cexCol = 0.45, cexRow = 0.45, main = "Sample-Sample Correlations",
          ColSideColors = pops) 

heatmapPlot

$rowInd
  [1] 342 341 340 339 338 337 336 335 334 333 332 331 330 329 328 327 326 325 324 323 322 321 320
 [24] 319 318 317 316 315 314 313 312 311 310 309 308 307 306 305 304 303 302 301 300 299 298 297
 [47] 296 295 294 293 292 291 290 289 288 287 286 285 284 283 282 281 280 279 278 277 276 275 274
 [70] 273 272 271 270 269 268 267 266 265 264 263 262 261 260 259 258 257 256 255 254 253 252 251
 [93] 250 249 248 247 246 245 244 243 242 241 240 239 238 237 236 235 234 233 232 231 230 229 228
[116] 227 226 225 224 223 222 221 220 219 218 217 216 215 214 213 212 211 210 209 208 207 206 205
[139] 204 203 202 201 200 199 198 197 196 195 194 193 192 191 190 189 188 187 186 185 184 183 182
[162] 181 180 179 178 177 176 175 174 173 172 171 170 169 168 167 166 165 164 163 162 161 160 159
[185] 158 157 156 155 154 153 152 151 150 149 148 147 146 145 144 143 142 141 140 139 138 137 136
[208] 135 134 133 132 131 130 129 128 127 126 125 124 123 122 121 120 119 118 117 116 115 114 113
[231] 112 111 110 109 108 107 106 105 104 103 102 101 100  99  98  97  96  95  94  93  92  91  90
[254]  89  88  87  86  85  84  83  82  81  80  79  78  77  76  75  74  73  72  71  70  69  68  67
[277]  66  65  64  63  62  61  60  59  58  57  56  55  54  53  52  51  50  49  48  47  46  45  44
[300]  43  42  41  40  39  38  37  36  35  34  33  32  31  30  29  28  27  26  25  24  23  22  21
[323]  20  19  18  17  16  15  14  13  12  11  10   9   8   7   6   5   4   3   2   1

$colInd
  [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23
 [24]  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46
 [47]  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69
 [70]  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92
 [93]  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
[116] 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
[139] 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
[162] 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
[185] 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
[208] 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
[231] 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
[254] 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
[277] 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
[300] 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322
[323] 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342

$call
heatmap.2(x = scaled_normalized_data_cor, Rowv = NA, Colv = NA, 
    dendrogram = "none", symm = T, col = mypal, trace = "none", 
    ColSideColors = pops, cexRow = 0.45, cexCol = 0.45, labRow = NA, 
    labCol = NA, main = "Sample-Sample Correlations")

$carpet
               HL60D7.M06_S89 HL60D7.M05_S80 HL60D7.M04_S69 HL60D7.M03_S61 HL60D7.M02_S21
HL60D0.C01_S1     -2.00000000    -2.00000000    -1.24622775    -1.15204907   -1.628853943
HL60D0.C02_S2     -2.00000000    -2.00000000    -1.23168502    -1.01412168   -1.644545859
               HL60D7.M01_S3 HL60D7.J02_S93 HL60D7.J01_S44 HL60D7.C87_S96 HL60D7.C86_S95
HL60D0.C01_S1    -1.05252571    -0.66865496   -1.712734386   -1.587968378   -0.521456275
HL60D0.C02_S2    -1.07088404    -0.89565369   -1.750494073   -2.000000000   -0.578798255
               HL60D7.C85_S94 HL60D7.C83_S91 HL60D7.C82_S90 HL60D7.C78_S85 HL60D7.C77_S84
HL60D0.C01_S1     -1.11530657    -1.71954101    0.470019132    -2.00000000   0.5631557388
HL60D0.C02_S2     -1.14160632    -1.68493727   -0.529114673    -2.00000000   0.3186543086
               HL60D7.C76_S83 HL60D7.C75_S82 HL60D7.C74_S81 HL60D7.C72_S78 HL60D7.C71_S77
HL60D0.C01_S1      0.38560090   -0.941742267   -0.287933412   0.0925599816   -0.966302866
HL60D0.C02_S2      0.49620243   -0.926914123   -0.335298192   0.0114647191   -0.910415156
               HL60D7.C70_S76 HL60D7.C69_S75 HL60D7.C68_S74 HL60D7.C66_S72 HL60D7.C65_S71
HL60D0.C01_S1     -1.28882932   0.2783873530    -0.21809895   0.5262184315    -1.64202515
HL60D0.C02_S2     -1.27114161   0.2208052703    -0.16814019   0.2361802986    -1.80718144
               HL60D7.C63_S68 HL60D7.C62_S67 HL60D7.C60_S64 HL60D7.C59_S63 HL60D7.C58_S62
HL60D0.C01_S1     -2.00000000   -1.489754697    -2.00000000    -1.08410341   -2.000000000
HL60D0.C02_S2     -2.00000000   -1.533454614    -2.00000000    -1.30825589   -2.000000000
               HL60D7.C56_S59 HL60D7.C55_S58 HL60D7.C53_S56 HL60D7.C52_S55 HL60D7.C51_S54
HL60D0.C01_S1     -1.40739065   -1.067748309   -1.814627411    -1.93942007    0.642718781
HL60D0.C02_S2     -1.49642640   -1.229317516   -1.885153426    -2.00000000    0.441304797
               HL60D7.C50_S53 HL60D7.C49_S52 HL60D7.C47_S50 HL60D7.C45_S48 HL60D7.C43_S46
HL60D0.C01_S1     -1.73972389    -0.70387149    -1.31236086   0.2097236057     -0.9810355
HL60D0.C02_S2     -1.73974620    -0.57706390    -1.33428723  -0.3888714340     -2.0000000
               HL60D7.C42_S45 HL60D7.C41_S43 HL60D7.C40_S42 HL60D7.C38_S40 HL60D7.C37_S39
HL60D0.C01_S1     -1.93718214    -1.57022512   -0.788605301    -1.39159290    -2.00000000
HL60D0.C02_S2     -2.00000000    -1.61970179   -0.688521879    -1.46255626    -2.00000000
               HL60D7.C36_S38 HL60D7.C34_S36 HL60D7.C33_S35 HL60D7.C32_S34 HL60D7.C31_S33
HL60D0.C01_S1     -0.82045931   -0.572937790    -1.48051343    0.079076352    -2.00000000
HL60D0.C02_S2     -0.84696042   -0.574296054    -1.54946038    0.106334236    -2.00000000
               HL60D7.C30_S32 HL60D7.C29_S31 HL60D7.C25_S27 HL60D7.C23_S25 HL60D7.C22_S24
HL60D0.C01_S1    -0.426397322    -0.57462361   -0.240715927   -0.023467820    -0.66437987
HL60D0.C02_S2    -0.206335894    -0.50608341   -0.401401774   -0.102051619    -0.58909278
               HL60D7.C20_S22 HL60D7.C19_S20 HL60D7.C18_S19 HL60D7.C15_S16 HL60D7.C14_S15
HL60D0.C01_S1   -0.5094498449   -1.815600977    -0.84106428   -1.521518353  -0.1908443722
HL60D0.C02_S2   -0.4804995794   -1.827682358    -0.82092114   -2.000000000  -0.1519123492
               HL60D7.C11_S12 HL60D7.C08_S9 HL60D7.C07_S8 HL60D7.C06_S7 HL60D7.C04_S5
HL60D0.C01_S1    4.497165e-01  -2.000000000  -0.783837001  -0.049563533    -1.8760975
HL60D0.C02_S2    1.946716e-01  -2.000000000  -0.717232887  -0.086377277    -1.9531983
               HL60D7.C03_S4 HL60D7.C01_S1 HL60D3.M05_S82 HL60D3.M04_S72 HL60D3.M03_S32
HL60D0.C01_S1  -0.0011016946   -0.84845137   -0.597062819  -0.9235468605    -1.58929221
HL60D0.C02_S2   0.0094505468   -1.22593853   -0.320485059  -0.6770440632    -1.49850971
               HL60D3.M02_S30 HL60D3.M01_S25 HL60D3.C89_S96 HL60D3.C88_S95 HL60D3.C87_S94
HL60D0.C01_S1     -1.96628613    -1.03595276    0.772162001  -0.7208413449   -0.495464374
HL60D0.C02_S2     -1.90205829    -0.87682871    0.494358188  -0.4876322206   -0.228507310
               HL60D3.C86_S93 HL60D3.C85_S92 HL60D3.C83_S90 HL60D3.C82_S89 HL60D3.C81_S88
HL60D0.C01_S1     -0.53347933    -1.50017564    -0.29688882   -0.624203384   -0.665133447
HL60D0.C02_S2     -0.31294890    -1.38753596    -0.13453576   -0.345152153   -0.313020648
               HL60D3.C80_S86 HL60D3.C79_S85 HL60D3.C78_S84 HL60D3.C77_S83 HL60D3.C76_S81
HL60D0.C01_S1    -0.436319163    -1.16314979    -0.88553322      0.5080606   -0.608856130
HL60D0.C02_S2    -0.177224469    -1.01570846    -0.71211681      0.7738598   -0.454030259
               HL60D3.C75_S80 HL60D3.C74_S79 HL60D3.C73_S78 HL60D3.C72_S77 HL60D3.C71_S76
HL60D0.C01_S1    -0.257084274    -1.26185054      0.4589765    -0.86357601   -0.681875090
HL60D0.C02_S2     0.131287251    -1.01279526      0.6252506    -0.63811150   -0.442899072
               HL60D3.C70_S75 HL60D3.C69_S74 HL60D3.C68_S73 HL60D3.C67_S71 HL60D3.C66_S70
HL60D0.C01_S1     0.179079469   -1.146197002   -1.367269638    -1.82432466    -0.39675766
HL60D0.C02_S2     0.336470950   -1.007353729   -1.256272212    -1.83352795    -0.18480254
               HL60D3.C65_S69 HL60D3.C64_S68 HL60D3.C63_S67 HL60D3.C62_S66 HL60D3.C61_S65
HL60D0.C01_S1    -2.000000000    -0.97164064  -0.4229960066    -1.46788004   -0.540434014
HL60D0.C02_S2    -2.000000000    -0.69232586  -0.1611924428    -1.48236971   -0.225230997
               HL60D3.C59_S63 HL60D3.C58_S62 HL60D3.C57_S61 HL60D3.C56_S60 HL60D3.C55_S59
HL60D0.C01_S1     0.152017592    -0.55481036   -0.606669745    -0.65223640   0.2557362006
HL60D0.C02_S2     0.324941870    -0.30281200   -0.344390005    -0.40708670   0.5305229375
               HL60D3.C54_S58 HL60D3.C53_S57 HL60D3.C52_S56 HL60D3.C51_S55 HL60D3.C50_S54
HL60D0.C01_S1    -0.910193847  -0.1062018447    -0.63667370   -0.653827176    -0.78479052
HL60D0.C02_S2    -0.800877520  -0.0246254614    -0.32020003   -0.343713232    -0.52478080
               HL60D3.C49_S53 HL60D3.C48_S52 HL60D3.C47_S51 HL60D3.C46_S50 HL60D3.C45_S49
HL60D0.C01_S1    -0.649432340    -1.39483673  -0.8916040567    -1.02238451    -0.52588783
HL60D0.C02_S2    -0.337301987    -1.22758932  -0.7252323072    -0.86722590    -0.36728241
               HL60D3.C44_S48 HL60D3.C43_S47 HL60D3.C42_S46 HL60D3.C41_S45 HL60D3.C40_S44
HL60D0.C01_S1     -1.07077852    -1.26867996     0.16768643    -1.02061931    -0.81840315
HL60D0.C02_S2     -0.96310961    -1.21234002     0.44393728    -0.80390751    -0.59363204
               HL60D3.C39_S43 HL60D3.C38_S42 HL60D3.C37_S41 HL60D3.C36_S40 HL60D3.C35_S39
HL60D0.C01_S1     -1.36437276   -0.505703621    -2.00000000    -1.28790196    -0.88494434
HL60D0.C02_S2     -1.23047416   -0.102275631    -2.00000000    -1.19350706    -0.67872153
               HL60D3.C34_S38 HL60D3.C32_S36 HL60D3.C31_S35 HL60D3.C30_S34 HL60D3.C29_S33
HL60D0.C01_S1     -1.14807843    -1.35444654    -1.85240914    -0.61256624    -1.34211482
HL60D0.C02_S2     -0.99710924    -1.19751812    -1.85777448    -0.49195114    -1.09473433
               HL60D3.C28_S31 HL60D3.C27_S29 HL60D3.C26_S28 HL60D3.C25_S27 HL60D3.C24_S26
HL60D0.C01_S1   -0.1113807749    -1.05580716   -1.092509362   -0.685542791    -1.13506108
HL60D0.C02_S2    0.2244249569    -0.81745135   -0.837334155   -0.460278942    -0.98230011
               HL60D3.C23_S24 HL60D3.C21_S22 HL60D3.C20_S21 HL60D3.C19_S20 HL60D3.C18_S19
HL60D0.C01_S1    -1.163903155    -1.70049703   -1.111428514     0.24928938    -2.00000000
HL60D0.C02_S2    -1.003617088    -1.74399062   -0.974161802     0.54942230    -2.00000000
               HL60D3.C17_S18 HL60D3.C16_S17 HL60D3.C15_S16 HL60D3.C14_S15 HL60D3.C13_S14
HL60D0.C01_S1     -1.01791262    -1.06302719   -0.581639135    -0.85951064   -1.348413638
HL60D0.C02_S2     -0.77250668    -0.80973589   -0.307520723    -0.67030214   -1.129733193
               HL60D3.C12_S13 HL60D3.C11_S12 HL60D3.C10_S11 HL60D3.C08_S9 HL60D3.C07_S8
HL60D0.C01_S1     -0.91373262    -1.94629718    -1.29027594   -1.52990325  -0.733605150
HL60D0.C02_S2     -0.76360287    -2.00000000    -1.12004258   -1.41675279  -0.475836010
               HL60D3.C06_S7 HL60D3.C05_S5 HL60D3.C04_S4 HL60D3.C03_S3 HL60D3.C02_S2
HL60D0.C01_S1    0.181229372  -1.584638530     0.5400507  -0.369017860  -0.165218871
HL60D0.C02_S2    0.377336090  -1.506579334     0.7703234  -0.152083379   0.172408632
               HL60D3.C01_S1 HL60D1.M01_S46 HL60D1.D02_S90 HL60D1.D01_S15 HL60D1.C92_S96
HL60D0.C01_S1   1.0018643796   0.3460524342    0.917237925   0.3977408549    1.124861185
HL60D0.C02_S2   0.7521304982   0.2602012325    0.391175320   0.3237837442    0.753310849
               HL60D1.C91_S95 HL60D1.C90_S94 HL60D1.C89_S93 HL60D1.C88_S92 HL60D1.C87_S91
HL60D0.C01_S1      0.90742945     0.96085546    0.841878470    0.887955444   0.8252383295
HL60D0.C02_S2      0.74194909     0.81755836    0.875057557    0.922747018   0.8968763793
               HL60D1.C85_S88 HL60D1.C84_S87 HL60D1.C83_S86 HL60D1.C82_S85 HL60D1.C81_S84
HL60D0.C01_S1     0.887819085     1.09886782     0.98023889    0.707563092    0.820051608
HL60D0.C02_S2     0.959006946     0.74807608     0.90777384    0.811796527    0.962380141
               HL60D1.C80_S83 HL60D1.C79_S82 HL60D1.C78_S81 HL60D1.C77_S80 HL60D1.C76_S79
HL60D0.C01_S1      0.80312442    1.041964638     0.87181730   0.5896751428     0.41232396
HL60D0.C02_S2      0.80160152    0.950545736     0.96243531   0.7518846578     0.67037433
               HL60D1.C75_S77 HL60D1.C74_S76 HL60D1.C73_S75 HL60D1.C72_S74 HL60D1.C71_S73
HL60D0.C01_S1     0.534541250    0.942340193    0.932425698   0.1503201019     0.19388997
HL60D0.C02_S2     0.673910529    0.920450177    0.798044621   0.2367959879     0.42919148
               HL60D1.C70_S72 HL60D1.C69_S71 HL60D1.C68_S70 HL60D1.C67_S69 HL60D1.C66_S68
HL60D0.C01_S1   -0.1345788304   0.8102082284   1.0006119860     0.93038771     0.95654370
HL60D0.C02_S2    0.1340436246   0.8213108335   0.8371979737     0.97888282     0.78971759
               HL60D1.C65_S67 HL60D1.C64_S66 HL60D1.C63_S65 HL60D1.C62_S64 HL60D1.C61_S63
HL60D0.C01_S1    1.0521559257     0.73464360    0.729014880     0.85333354     1.00218452
HL60D0.C02_S2    0.8649270273     0.74490752    0.887515387     0.84079398     0.84442060
               HL60D1.C60_S62 HL60D1.C59_S61 HL60D1.C58_S60 HL60D1.C57_S59 HL60D1.C56_S58
HL60D0.C01_S1      0.60515940    0.006853455    0.404461062     0.89504593     0.82194874
HL60D0.C02_S2      0.79708381   -0.550320796    0.645550675     0.88167915     0.78483473
               HL60D1.C55_S57 HL60D1.C54_S56 HL60D1.C53_S55 HL60D1.C52_S54 HL60D1.C51_S53
HL60D0.C01_S1    0.0956060000    0.912045073     0.90374242     0.96022661     0.80948764
HL60D0.C02_S2    0.1491292769    0.863215479     0.91915593     0.90584245     0.85129974
               HL60D1.C50_S52 HL60D1.C49_S51 HL60D1.C48_S50 HL60D1.C47_S49 HL60D1.C46_S48
HL60D0.C01_S1     0.880597248    0.780644844     0.77993131    0.404975170    0.506026227
HL60D0.C02_S2     0.888912417    0.803352265     0.82966400    0.572735645    0.518209115
               HL60D1.C45_S47 HL60D1.C44_S45 HL60D1.C43_S44 HL60D1.C42_S43 HL60D1.C41_S42
HL60D0.C01_S1      0.88973692    0.895607468    0.774632036     0.94683416   0.7605866731
HL60D0.C02_S2      0.89287166    0.976607355    0.925026730     0.75194174   0.8626161189
               HL60D1.C40_S41 HL60D1.C39_S40 HL60D1.C38_S39 HL60D1.C37_S38 HL60D1.C36_S37
HL60D0.C01_S1       0.7112114   0.9108045402    0.809994052   0.6894151076    0.848631651
HL60D0.C02_S2       0.7565713   0.9954046814    0.366112051   0.7936724452    0.957635252
               HL60D1.C34_S35 HL60D1.C33_S34 HL60D1.C32_S33 HL60D1.C31_S32 HL60D1.C30_S31
HL60D0.C01_S1     0.688935003     0.94883087    0.904772370    0.310670729     0.65721093
HL60D0.C02_S2     0.184711143     0.73595236    0.954423748    0.570880286     0.17588279
               HL60D1.C29_S30 HL60D1.C28_S29 HL60D1.C27_S28 HL60D1.C26_S27 HL60D1.C25_S26
HL60D0.C01_S1     0.800206370     0.73387061     0.90349578    0.799510224     0.90228084
HL60D0.C02_S2     0.846652838     0.80086184     0.95895959    0.918590274     0.64615833
               HL60D1.C24_S25 HL60D1.C23_S24 HL60D1.C22_S23 HL60D1.C21_S22 HL60D1.C20_S21
HL60D0.C01_S1    0.3522433265    0.828483929     1.03460555     0.93586209   -0.316419861
HL60D0.C02_S2    0.5782122081    0.952791016     0.90975092     0.71497730   -0.249072450
               HL60D1.C19_S20 HL60D1.C18_S19 HL60D1.C17_S18 HL60D1.C16_S17 HL60D1.C15_S16
HL60D0.C01_S1      0.97859494     0.94202926    0.837149994     0.76210439    0.549792998
HL60D0.C02_S2      1.01721889     0.75266736    0.815890779     0.94583979    0.665877229
               HL60D1.C14_S14 HL60D1.C13_S13 HL60D1.C12_S12 HL60D1.C11_S11 HL60D1.C10_S10
HL60D0.C01_S1      0.96311644     1.03288232    0.909852378    0.975617065    0.912639846
HL60D0.C02_S2      1.00218787     0.87312988    0.897382361    0.989007075    0.718396555
               HL60D1.C09_S9 HL60D1.C08_S8 HL60D1.C07_S7 HL60D1.C06_S6 HL60D1.C05_S5
HL60D0.C01_S1     0.90326258   0.745645089    0.58861939    0.88056616   0.853567429
HL60D0.C02_S2     0.79571790   0.907236005   -0.04630890    0.91728855   0.893584614
               HL60D1.C04_S4 HL60D1.C03_S3 HL60D1.C02_S2 HL60D1.C01_S1 HL60D0.M01_S40
HL60D0.C01_S1    0.945883486  0.8900318441     0.4004916  1.1050638538   -0.023218635
HL60D0.C02_S2    0.935603286  0.9404398426     0.6524690  0.4887278258    0.331577487
               HL60D0.C94_S96 HL60D0.C93_S95 HL60D0.C92_S94 HL60D0.C91_S93 HL60D0.C90_S92
HL60D0.C01_S1    1.1394201062     0.82951278     0.89298688    1.049068198   0.5261074573
HL60D0.C02_S2    0.8108235067     0.70653961     0.56809681    0.734517656   0.8770376334
               HL60D0.C89_S91 HL60D0.C88_S90 HL60D0.C87_S89 HL60D0.C86_S88 HL60D0.C85_S87
HL60D0.C01_S1      0.88359949     0.83525075     0.80978696    0.602256954     0.90519534
HL60D0.C02_S2      0.93190699     0.27160573     0.85856481    0.834108053     0.26874235
               HL60D0.C84_S86 HL60D0.C83_S85 HL60D0.C82_S84 HL60D0.C81_S83 HL60D0.C80_S82
HL60D0.C01_S1     0.888867288     0.83846517   0.6925914817     0.92447498    0.875684423
HL60D0.C02_S2     0.943415035     0.38979991   0.6667079288     0.76733055    0.986477354
               HL60D0.C79_S81 HL60D0.C78_S80 HL60D0.C77_S79 HL60D0.C76_S78 HL60D0.C75_S77
HL60D0.C01_S1     0.940563952     0.91007169    0.579296815     0.91094050    0.731764464
HL60D0.C02_S2     0.945640995     0.97018763    0.469078060     0.49419144    0.306801128
               HL60D0.C74_S76 HL60D0.C73_S75 HL60D0.C72_S74 HL60D0.C71_S73 HL60D0.C69_S71
HL60D0.C01_S1      0.90099421     1.05983983     0.81575544    0.944878951    0.846767882
HL60D0.C02_S2      0.81666394     0.79880121     0.56856759    1.033410660    0.468232683
               HL60D0.C68_S70 HL60D0.C67_S69 HL60D0.C65_S67 HL60D0.C64_S66 HL60D0.C63_S65
HL60D0.C01_S1    -0.276847711      0.5240458    0.888430922     0.98299231    0.630523954
HL60D0.C02_S2    -0.033350076      0.6996085    0.818349412     1.00170439    0.380655740
               HL60D0.C62_S64 HL60D0.C61_S63 HL60D0.C60_S62 HL60D0.C59_S61 HL60D0.C58_S60
HL60D0.C01_S1      0.86861404    0.930233377     0.60693934    0.884353985     0.77138766
HL60D0.C02_S2      0.86302797    0.972857816     0.48148709    1.028962115     0.47273024
               HL60D0.C57_S59 HL60D0.C56_S58 HL60D0.C55_S57 HL60D0.C54_S56 HL60D0.C53_S55
HL60D0.C01_S1    0.2513210003    0.914418028    0.832649923    0.562595703    0.930483678
HL60D0.C02_S2    0.4563382556    0.964086191    0.991154331    0.814203145    0.999747427
               HL60D0.C52_S54 HL60D0.C51_S53 HL60D0.C50_S52 HL60D0.C49_S51 HL60D0.C48_S50
HL60D0.C01_S1     0.762608395   0.7949086898    0.742013817    0.907624301   0.7844687097
HL60D0.C02_S2     0.904173917   0.8131753318    0.548211323    0.866015882   1.0078908745
               HL60D0.C47_S49 HL60D0.C46_S48 HL60D0.C45_S47 HL60D0.C44_S46 HL60D0.C43_S45
HL60D0.C01_S1      0.94674615    0.743026687    0.464038894    0.488063469     0.69929306
HL60D0.C02_S2      1.01387187    0.948777050    0.679185520    0.472655287     0.90748332
               HL60D0.C42_S44 HL60D0.C41_S43 HL60D0.C40_S42 HL60D0.C39_S41 HL60D0.C38_S39
HL60D0.C01_S1      0.84833311   0.8383853954    0.708591878    -1.16543188    0.811806304
HL60D0.C02_S2      0.92963880   0.9645451695    0.761574784    -1.03023555    0.934231165
               HL60D0.C37_S38 HL60D0.C36_S37 HL60D0.C35_S36 HL60D0.C34_S35 HL60D0.C33_S34
HL60D0.C01_S1      0.53797356   0.6869112713     0.66998421     0.74866574    0.941145188
HL60D0.C02_S2      0.43094313   0.4172068952     0.69860524     0.60987486    1.048529989
               HL60D0.C32_S33 HL60D0.C31_S32 HL60D0.C30_S31 HL60D0.C29_S30 HL60D0.C28_S29
HL60D0.C01_S1      0.50676735    0.409480981    0.977930000     0.60544852     0.97577993
HL60D0.C02_S2      0.52855249    0.703768578    0.955693595     0.86197999     0.93367492
               HL60D0.C27_S28 HL60D0.C26_S27 HL60D0.C25_S26 HL60D0.C24_S25 HL60D0.C23_S24
HL60D0.C01_S1     0.815346297    0.889187245    0.831616917    0.908100454   0.9225887397
HL60D0.C02_S2     0.928519674    0.680932740    0.634464647    0.891497921   1.0042573513
               HL60D0.C22_S23 HL60D0.C21_S22 HL60D0.C20_S21 HL60D0.C19_S20 HL60D0.C18_S19
HL60D0.C01_S1     0.727526318     0.59252512     0.81817460    0.878279878     0.92022847
HL60D0.C02_S2     0.921481309     0.66186244     0.62392661    1.009671626     0.48073024
               HL60D0.C17_S18 HL60D0.C16_S17 HL60D0.C15_S16 HL60D0.C14_S15 HL60D0.C13_S14
HL60D0.C01_S1      0.95063271    0.778792926     0.94684080    0.219167326   7.417246e-01
HL60D0.C02_S2      0.91159137    0.957939265     0.99849159    0.271894550   5.392685e-01
               HL60D0.C12_S13 HL60D0.C11_S12 HL60D0.C10_S11 HL60D0.C09_S10 HL60D0.C08_S8
HL60D0.C01_S1      0.85977009    0.358939583    0.917812986     1.04561781   0.461830655
HL60D0.C02_S2      0.60883213    0.668507078    0.981373981     0.56827850   0.220483784
               HL60D0.C07_S7 HL60D0.C06_S6 HL60D0.C05_S5 HL60D0.C04_S4 HL60D0.C03_S3
HL60D0.C01_S1    0.722335912    0.89777844    0.86303856   0.711937944    0.88209417
HL60D0.C02_S2    0.004756029    0.66771148    0.97401816  -0.105555220    0.74689043
               HL60D0.C02_S2 HL60D0.C01_S1
HL60D0.C01_S1    0.967421282   1.202032179
HL60D0.C02_S2    1.091237000   0.771333522
 [ reached getOption("max.print") -- omitted 340 rows ]

$rowDendrogram
'dendrogram' with 2 branches and 342 members total, at height 1.414214 

$colDendrogram
'dendrogram' with 2 branches and 342 members total, at height 1.414214 

$breaks
  [1] -2.00 -1.96 -1.92 -1.88 -1.84 -1.80 -1.76 -1.72 -1.68 -1.64 -1.60 -1.56 -1.52 -1.48 -1.44
 [16] -1.40 -1.36 -1.32 -1.28 -1.24 -1.20 -1.16 -1.12 -1.08 -1.04 -1.00 -0.96 -0.92 -0.88 -0.84
 [31] -0.80 -0.76 -0.72 -0.68 -0.64 -0.60 -0.56 -0.52 -0.48 -0.44 -0.40 -0.36 -0.32 -0.28 -0.24
 [46] -0.20 -0.16 -0.12 -0.08 -0.04  0.00  0.04  0.08  0.12  0.16  0.20  0.24  0.28  0.32  0.36
 [61]  0.40  0.44  0.48  0.52  0.56  0.60  0.64  0.68  0.72  0.76  0.80  0.84  0.88  0.92  0.96
 [76]  1.00  1.04  1.08  1.12  1.16  1.20  1.24  1.28  1.32  1.36  1.40  1.44  1.48  1.52  1.56
 [91]  1.60  1.64  1.68  1.72  1.76  1.80  1.84  1.88  1.92  1.96  2.00

$col
  [1] "#053061" "#073568" "#0A3A70" "#0D4077" "#10457F" "#134B86" "#15508E" "#185696" "#1B5B9D"
 [10] "#1E61A5" "#2166AC" "#246AAE" "#286FB0" "#2B74B3" "#2F78B5" "#327DB7" "#3581BA" "#3986BC"
 [19] "#3C8ABE" "#408FC1" "#4494C3" "#4C99C6" "#549EC9" "#5CA3CB" "#64A8CE" "#6CADD1" "#74B2D3"
 [28] "#7CB7D6" "#84BCD9" "#8CC1DC" "#93C5DE" "#9AC9E0" "#A0CCE2" "#A7CFE4" "#ADD2E5" "#B3D6E7"
 [37] "#BAD9E9" "#C0DCEB" "#C6DFED" "#CDE3EE" "#D2E5F0" "#D6E7F0" "#DAE9F1" "#DEEBF2" "#E1EDF3"
 [46] "#E5EEF3" "#E9F0F4" "#EDF2F5" "#F1F4F5" "#F5F6F6" "#F7F5F4" "#F7F2EF" "#F8EFEA" "#F9EDE6"
 [55] "#F9EAE1" "#FAE7DC" "#FAE4D7" "#FBE1D2" "#FCDECD" "#FCDCC8" "#FCD7C2" "#FBD2BB" "#FACCB4"
 [64] "#F9C7AD" "#F8C1A6" "#F7BC9F" "#F7B799" "#F6B192" "#F5AC8B" "#F4A684" "#F1A07E" "#EE9978"
 [73] "#EB9273" "#E88B6E" "#E58468" "#E27D63" "#DF765E" "#DC6F58" "#D96853" "#D6614E" "#D35A4A"
 [82] "#CF5246" "#CB4B43" "#C8443F" "#C43D3C" "#C03539" "#BD2E35" "#B92732" "#B6202E" "#B2182B"
 [91] "#AB1529" "#A31328" "#9C1027" "#940E26" "#8C0C25" "#850923" "#7D0722" "#760421" "#6E0220"
[100] "#67001F"

$colorTable

$layout
$layout$lmat
     [,1] [,2]
[1,]    5    4
[2,]    0    1
[3,]    3    2

$layout$lhei
[1] 1.5 0.2 4.0

$layout$lwid
[1] 1.5 4.0

#will not be needed: don't keep big data
remove(normalized_data_cor)
remove(scaled_normalized_data_cor)
#remove(heatmapPlot)

In this heat map we clearly see the high correlation between the samples of day 0 and 1 and between the samples of day 3 and 7.

t-SNE plots

General t-SNE

Let’s plot the t-SNE plots of the cells that aren’t dead to see if we can find clusters and / or outliers. Note that we aren’t taking into account genes that aren’t expressed in any samples.

# t-SNE plot with alive cells and non zero genes
#ordering data by expn
alive_tsne_data <- cbind(dataNZ, "avg_expn" = rowMeans(dataNZ), "miRNA" = row.names(dataNZ)) 
alive_tsne_data <- dplyr::arrange(alive_tsne_data,desc(avg_expn))
alive_tsne_data <- t(alive_tsne_data)
colnames(alive_tsne_data) <- alive_tsne_data["miRNA", ]
#removes rows that were added for ordering
alive_tsne_data <- alive_tsne_data[!rownames(alive_tsne_data) %in% c("miRNA","avg_expn"), ]
tsne_out <- Rtsne(dist(alive_tsne_data)) # Run TSNE
# Show the objects in the 2D tsne representation
tSNETotal <- qplot(x=tsne_out$Y[,1], y=tsne_out$Y[,2],color=data_day) + 
  labs(colour = "Cell Type", x = "tsne1", y = "tsne2")+
  ggtitle("t-SNE alive cells and non zero genes") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
tSNETotal

#will not be needed: don't keep big data
remove(tsne_out)
#remove(tSNETotal)

As with the correlation heatmap, we see that the samples from day 0 and 1 are clustered together while those from day 3 and 7 are clustered together.

t-SNE of Day 0 and Day 1

Now that we have seen that the samples from day 0 and 1 are closely related. Let’s try to take only those 2 in order to see if we can separate them from each other.

#tsne on Day 0/Day 1 cells only
indexDay01 <- which(data_day %in% c("Day 0","Day 1"))
group_one_data <- alive_tsne_data[indexDay01, ] #selecting only day 0 and day 1 cells
tsne_out <- Rtsne(dist(group_one_data)) # Run TSNE
# Show the objects in the 2D tsne representation
tsneDay01 <- qplot(x=tsne_out$Y[,1], y=tsne_out$Y[,2],color=data_day[indexDay01]) +
  labs(colour = "Cell Type", x = "tsne1", y = "tsne2") +
  ggtitle("t-SNE only day 0 and 1") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
tsneDay01

#will not be needed: don't keep big data
remove(tsne_out)
#remove(tsneDay01)

The t-SNE can clearly not separate the 2 samples from each other. Finally let’s investigate whether they can be separated with a t-SNE in 3 dimension:

tsne_out <- Rtsne(dist(group_one_data), dims = 3) # Run TSNE in 3 dimensions
plot_ly(x=tsne_out$Y[,1],
        y=tsne_out$Y[,2],
        z=tsne_out$Y[,3],
        type="scatter3d",
        color=data_day[indexDay01],
        mode="markers") %>% 
  layout(title = 't-SNE day 0 and 1 3D')

#will not be needed: don't keep big data
remove(tsne_out)
remove(group_one_data)
#remove(tsneDay01)

We can still not distinguish both type of cells. This really shows how similar they are.

t-SNE of Day 0 and Day 1

Now let’s look at the cluster Day 3 and 7, to see if they can be distinguished.

#tsne on Day 0/Day 1 cells only
indexDay37 <- which(data_day %in% c("Day 3","Day 7"))
group_two_data <- alive_tsne_data[indexDay37, ] #selecting only day 0 and day 1 cells
tsne_out <- Rtsne(dist(group_two_data)) # Run TSNE
# Show the objects in the 2D tsne representation
tsneDay37 <- qplot(x=tsne_out$Y[,1], y=tsne_out$Y[,2],color=data_day[indexDay37]) +
  labs(colour = "Cell Type", x = "tsne1", y = "tsne2") +
  ggtitle("t-SNE only day 3 and 7") +
  theme(text = element_text(size=14),plot.title = element_text(hjust = 0.5)) 
tsneDay37

#will not be needed: don't keep big data
remove(tsne_out)
#remove(tsneDay37)

As before, the 2 type of cells cannot be distinguished with a t-SNE plot in 2D. They will probably not be distinguishable in 3D but let’s be sure:

tsne_out <- Rtsne(dist(group_two_data), dims = 3) # Run TSNE in 3 dimensions
plot_ly(x=tsne_out$Y[,1],
        y=tsne_out$Y[,2],
        z=tsne_out$Y[,3],
        type="scatter3d",
        color=data_day[indexDay37],
        mode="markers") %>% 
  layout(title = 't-SNE day 3 and 7 3D')

#will not be needed: don't keep big data
remove(tsne_out)
remove(group_two_data)

Not suprisingly, the second cluster cannot be seprated further with 3 dimensional t-SNE.

Fit linear model

Note that for fitting linear models, we cannot use the DATASEQ library as it requires counts (i.e. discrete values). We now have continuous values because of they way we normalized the data.

Voom

Mean-Variance trend

First let us visualize the mean-variance trend, to see. Note: here it is important that we use only the genes with no non zero values for every sample as we can see in paragraph 15.3 of (limma manual)[https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf]. From the same article they say that as a rule of thumb we should keep only genes that have a

#working ... HOW TO USE DATADESEQ ON NORMALIZED => NOT COUNTS ANY MORE?? + IS IT WORTH IT ? cAN we round ?
#dataDESeq <- rlog(as.matrix(data), blind=FALSE )
#makes a convenient dge object for analysis
dge <- DGEList(counts=dataNZ, group = data_day)
# applies TMM normalization to dge
#isExpr <- rowSums(dge$counts) > 1000
#dge <- dge[isExpr,]
dge <- calcNormFactors(dge)
data_voomed <- voom(dge,design,plot=TRUE)

#will not be needed: don't keep big data
remove(dge)

We wee that the mean variance trend seam to be decreasing very quickly and steadily, we can thus conclude that the experiment seem to be of low biological variation. For more information please see the “Removing heteroscedascity from count data” subpart of (this paper)[https://f1000research.com/articles/5-1408]

Here is an example of the data we have:

data_voomed[1:2,1:10]

An object of class "EList"
$targets
               group  lib.size norm.factors
HL60D0.C01_S1  Day 0 100862.30    0.9833482
HL60D0.C02_S2  Day 0  68768.90    0.9420345
HL60D0.C03_S3  Day 0  60552.35    1.0147531
HL60D0.C04_S4  Day 0 279868.19    0.9569073
HL60D0.C05_S5  Day 0  65361.77    1.1399012
HL60D0.C06_S6  Day 0  36823.33    0.8372045
HL60D0.C07_S7  Day 0 127944.17    0.8744118
HL60D0.C08_S8  Day 0  51630.86    0.9873263
HL60D0.C09_S10 Day 0 217824.43    1.0843357
HL60D0.C10_S11 Day 0  55444.19    1.1485957

$E
              HL60D0.C01_S1 HL60D0.C02_S2 HL60D0.C03_S3 HL60D0.C04_S4 HL60D0.C05_S5 HL60D0.C06_S6
hsa-let-7a-5p     13.374469     12.879195      12.91583     7.6516046      13.06635      13.48863
hsa-let-7a-3p      2.309527      2.862079      10.11151     0.8371754      10.14368      10.50555
              HL60D0.C07_S7 HL60D0.C08_S8 HL60D0.C09_S10 HL60D0.C10_S11
hsa-let-7a-5p     13.102009      13.73708       12.72418      13.310086
hsa-let-7a-3p      1.966402      10.01462        9.35818       9.870383

$weights
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]      [,9]
[1,] 0.5885646 0.5004936 0.4749981 0.9426669 0.4901102 0.3899081 0.6533318 0.4453390 0.8348369
[2,] 0.2090121 0.1851257 0.1774882 0.2881312 0.1820762 0.1454568 0.2249386 0.1678079 0.2654087
         [,10]
[1,] 0.4582897
[2,] 0.1720934

$design
   HL60D0 HL60D1 HL60D3 HL60D7
1       1      0      0      0
2       1      0      0      0
3       1      0      0      0
4       1      0      0      0
5       1      0      0      0
6       1      0      0      0
7       1      0      0      0
8       1      0      0      0
9       1      0      0      0
10      1      0      0      0

Fit linear model

fit <- lmFit(data_voomed, design)
fit <- contrasts.fit(fit, contrasts=contrastMatrix)
efit <- eBayes(fit)
topTable(efit, coef=colnames(coef(efit)))

plotSA(efit)

#will not be needed: don't keep big data
remove(fit)
remove(data_voomed)

We can thus see from this plot that the variance is no longer dependent on the mean expression level.

Examining the number of DE genes

Let’s now look at differential expression levels of genes. Note that the adjusted p-value cutoff is 5%.

dt <- decideTests(efit)
summary(dt)

   HL60D1 - HL60D0 HL60D3 - HL60D0 HL60D7 - HL60D0 HL60D1 - HL60D3 HL60D1 - HL60D7
-1              28             325             498              87              41
0              705             376             206             292             221
1               11              43              40             365             482
   HL60D3 - HL60D7
-1              24
0              624
1               96

As we have already previously seen in the multiple plots of the previous section: most differential expressed genes are between day 0/1 and 3/7. We also note that a lot of genes in Day 3 and 7 are downregulated (compared to day 0 and 1).

Finding genes that are often differentially expressed compared to day 0

Let’s look at the genes that are always differentially expressed compared to day 0

# Finding genes that are not zero in every comparaison with day 0
commonDEgenesDay0 <- which(dt[,1]!=0 & dt[,2]!=0 & dt[,3]!=0)
names(commonDEgenesDay0)

 [1] "hsa-let-7a-3p"   "hsa-let-7b-3p"   "hsa-miR-1260a"   "hsa-miR-132-3p"  "hsa-miR-132-5p" 
 [6] "hsa-miR-143-3p"  "hsa-miR-17-3p"   "hsa-miR-181a-3p" "hsa-miR-181b-5p" "hsa-miR-181b-3p"
[11] "hsa-miR-199b-5p" "hsa-miR-19b-3p"  "hsa-miR-28-3p"   "hsa-miR-30d-5p"  "hsa-miR-5701"   
[16] "hsa-miR-7974"   

Let’s visualise it with a nice venn diagramm

vennDiagram(dt[,1:3], circle.col=c("turquoise", "salmon", "green"))

Finding genes that are often differentially between day 0/1 and day 3/7

Let’s look at the genes that are always differentially between day 0/1 and 3/7 as they seem to be the ones of interest.

# Finding genes that are not zero in every comparaison with day 0
commonDEgenesDay0137 <- which(dt[,2]!=0 & dt[,3]!=0 & dt[,4]!=0 & dt[,5]!=0)
names(commonDEgenesDay0137)

  [1] "hsa-let-7c-3p"     "hsa-let-7d-3p"     "hsa-let-7e-5p"     "hsa-let-7i-5p"    
  [5] "hsa-miR-101-3p"    "hsa-miR-106b-3p"   "hsa-miR-106b-5p"   "hsa-miR-1233-3p"  
  [9] "hsa-miR-1234-3p"   "hsa-miR-1236-3p"   "hsa-miR-125a-5p"   "hsa-miR-125b-2-3p"
 [13] "hsa-miR-1260a"     "hsa-miR-1268a"     "hsa-miR-1268b"     "hsa-miR-1271-5p"  
 [17] "hsa-miR-1271-3p"   "hsa-miR-1273d"     "hsa-miR-1273h-5p"  "hsa-miR-1275"     
 [21] "hsa-miR-1278"      "hsa-miR-1281"      "hsa-miR-1285-3p"   "hsa-miR-1287-5p"  
 [25] "hsa-miR-129-2-3p"  "hsa-miR-1291"      "hsa-miR-1303"      "hsa-miR-1304-3p"  
 [29] "hsa-miR-1306-3p"   "hsa-miR-1307-3p"   "hsa-miR-132-3p"    "hsa-miR-138-5p"   
 [33] "hsa-miR-143-5p"    "hsa-miR-143-3p"    "hsa-miR-145-5p"    "hsa-miR-146b-3p"  
 [37] "hsa-miR-149-5p"    "hsa-miR-150-5p"    "hsa-miR-151a-5p"   "hsa-miR-152-5p"   
 [41] "hsa-miR-1537-5p"   "hsa-miR-155-3p"    "hsa-miR-181a-3p"   "hsa-miR-181b-3p"  
 [45] "hsa-miR-181c-3p"   "hsa-miR-182-5p"    "hsa-miR-1827"      "hsa-miR-186-3p"   
 [49] "hsa-miR-190a-5p"   "hsa-miR-192-5p"    "hsa-miR-193a-5p"   "hsa-miR-193a-3p"  
 [53] "hsa-miR-195-5p"    "hsa-miR-196a-5p"   "hsa-miR-196b-5p"   "hsa-miR-197-5p"   
 [57] "hsa-miR-19b-3p"    "hsa-miR-200c-5p"   "hsa-miR-205-5p"    "hsa-miR-210-3p"   
 [61] "hsa-miR-2115-3p"   "hsa-miR-2116-3p"   "hsa-miR-212-3p"    "hsa-miR-218-5p"   
 [65] "hsa-miR-223-3p"    "hsa-miR-224-5p"    "hsa-miR-2277-3p"   "hsa-miR-2277-5p"  
 [69] "hsa-miR-2278"      "hsa-miR-23a-5p"    "hsa-miR-23a-3p"    "hsa-miR-25-3p"    
 [73] "hsa-miR-25-5p"     "hsa-miR-27a-5p"    "hsa-miR-29b-1-5p"  "hsa-miR-301b-5p"  
 [77] "hsa-miR-3074-3p"   "hsa-miR-3115"      "hsa-miR-3128"      "hsa-miR-3136-5p"  
 [81] "hsa-miR-3146"      "hsa-miR-3154"      "hsa-miR-3158-3p"   "hsa-miR-3173-3p"  
 [85] "hsa-miR-3177-3p"   "hsa-miR-3180-5p"   "hsa-miR-3181"      "hsa-miR-3198"     
 [89] "hsa-miR-328-3p"    "hsa-miR-338-5p"    "hsa-miR-33b-3p"    "hsa-miR-34a-5p"   
 [93] "hsa-miR-3529-3p"   "hsa-miR-3607-5p"   "hsa-miR-3607-3p"   "hsa-miR-3614-3p"  
 [97] "hsa-miR-3648"      "hsa-miR-365a-3p"   "hsa-miR-365b-3p"   "hsa-miR-3662"     
[101] "hsa-miR-3677-3p"   "hsa-miR-3680-5p"   "hsa-miR-3680-3p"   "hsa-miR-3691-5p"  
[105] "hsa-miR-378b"      "hsa-miR-378g"      "hsa-miR-3912-5p"   "hsa-miR-3913-5p"  
[109] "hsa-miR-3917"      "hsa-miR-3922-3p"   "hsa-miR-3940-3p"   "hsa-miR-3940-5p"  
[113] "hsa-miR-3942-5p"   "hsa-miR-3944-3p"   "hsa-miR-409-3p"    "hsa-miR-422a"     
[117] "hsa-miR-424-3p"    "hsa-miR-425-5p"    "hsa-miR-4284"      "hsa-miR-4421"     
[121] "hsa-miR-4426"      "hsa-miR-4429"      "hsa-miR-4434"      "hsa-miR-4435"     
[125] "hsa-miR-4449"      "hsa-miR-4454"      "hsa-miR-4477b"     "hsa-miR-4484"     
[129] "hsa-miR-450a-2-3p" "hsa-miR-4510"      "hsa-miR-4515"      "hsa-miR-4518"     
[133] "hsa-miR-4523"      "hsa-miR-455-5p"    "hsa-miR-4660"      "hsa-miR-4662a-5p" 
[137] "hsa-miR-4664-5p"   "hsa-miR-4685-3p"   "hsa-miR-4707-5p"   "hsa-miR-4723-5p"  
[141] "hsa-miR-4743-3p"   "hsa-miR-4744"      "hsa-miR-4745-5p"   "hsa-miR-4745-3p"  
[145] "hsa-miR-4747-3p"   "hsa-miR-4766-3p"   "hsa-miR-4772-3p"   "hsa-miR-4777-3p"  
[149] "hsa-miR-4779"      "hsa-miR-4782-3p"   "hsa-miR-4785"      "hsa-miR-4791"     
[153] "hsa-miR-4796-5p"   "hsa-miR-4797-3p"   "hsa-miR-4802-3p"   "hsa-miR-486-5p"   
[157] "hsa-miR-495-3p"    "hsa-miR-497-5p"    "hsa-miR-5001-5p"   "hsa-miR-5004-3p"  
[161] "hsa-miR-5094"      "hsa-miR-5188"      "hsa-miR-548a-3p"   "hsa-miR-548aa"    
[165] "hsa-miR-548ag"     "hsa-miR-548am-5p"  "hsa-miR-548ar-3p"  "hsa-miR-548au-3p" 
[169] "hsa-miR-548aw"     "hsa-miR-548ay-5p"  "hsa-miR-548b-5p"   "hsa-miR-548d-5p"  
[173] "hsa-miR-548f-3p"   "hsa-miR-548n"      "hsa-miR-548s"      "hsa-miR-548t-5p"  
[177] "hsa-miR-548t-3p"   "hsa-miR-550b-2-5p" "hsa-miR-552-5p"    "hsa-miR-5585-5p"  
[181] "hsa-miR-5586-5p"   "hsa-miR-5684"      "hsa-miR-5690"      "hsa-miR-570-3p"   
[185] "hsa-miR-573"       "hsa-miR-581"       "hsa-miR-582-5p"    "hsa-miR-597-3p"   
[189] "hsa-miR-6134"      "hsa-miR-617"       "hsa-miR-628-3p"    "hsa-miR-629-5p"   
[193] "hsa-miR-639"       "hsa-miR-6500-3p"   "hsa-miR-6509-5p"   "hsa-miR-6513-3p"  
[197] "hsa-miR-663a"      "hsa-miR-663b"      "hsa-miR-6724-5p"   "hsa-miR-6726-3p"  
[201] "hsa-miR-6727-5p"   "hsa-miR-6729-3p"   "hsa-miR-6730-5p"   "hsa-miR-6733-3p"  
[205] "hsa-miR-6734-3p"   "hsa-miR-6734-5p"   "hsa-miR-6735-3p"   "hsa-miR-6737-3p"  
[209] "hsa-miR-6749-3p"   "hsa-miR-6753-5p"   "hsa-miR-6754-3p"   "hsa-miR-6762-3p"  
[213] "hsa-miR-6763-5p"   "hsa-miR-6764-3p"   "hsa-miR-6764-5p"   "hsa-miR-6774-3p"  
[217] "hsa-miR-6777-3p"   "hsa-miR-6781-3p"   "hsa-miR-6781-5p"   "hsa-miR-6787-3p"  
[221] "hsa-miR-6790-3p"   "hsa-miR-6793-3p"   "hsa-miR-6797-3p"   "hsa-miR-6798-3p"  
[225] "hsa-miR-6803-3p"   "hsa-miR-6812-3p"   "hsa-miR-6813-5p"   "hsa-miR-6818-3p"  
[229] "hsa-miR-6819-3p"   "hsa-miR-6833-3p"   "hsa-miR-6837-5p"   "hsa-miR-6844"     
[233] "hsa-miR-6849-3p"   "hsa-miR-6852-5p"   "hsa-miR-6857-3p"   "hsa-miR-6858-5p"  
[237] "hsa-miR-6859-3p"   "hsa-miR-6862-3p"   "hsa-miR-6868-3p"   "hsa-miR-6873-3p"  
[241] "hsa-miR-6885-3p"   "hsa-miR-6889-3p"   "hsa-miR-6891-5p"   "hsa-miR-6891-3p"  
[245] "hsa-miR-6894-3p"   "hsa-miR-7111-3p"   "hsa-miR-7112-3p"   "hsa-miR-7114-3p"  
[249] "hsa-miR-7641"      "hsa-miR-766-5p"    "hsa-miR-7705"      "hsa-miR-7854-3p"  
[253] "hsa-miR-8485"      "hsa-miR-874-3p"    "hsa-miR-877-5p"    "hsa-miR-885-5p"   
[257] "hsa-miR-9-3p"      "hsa-miR-92a-1-5p"  "hsa-miR-92a-3p"    "hsa-miR-937-5p"   
[261] "hsa-miR-940"       "hsa-miR-942-3p"    "hsa-miR-99a-3p"   

Let’s visualise it with a nice venn diagramm

vennDiagram(dt[,2:5], circle.col=c("turquoise", "salmon", "green","yellow"))